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Abstract: 

This paper focuses on investigation of millimeter wave (mmW) extinction due to its interaction with 2-D 
layers of randomly distributed air-bubbles on the surface of the ocean. The Split-step Fourier method was 
adopted for evaluation of the refractive and diffractive effects of scattered mmW due to its interaction with 
layers of foam-covered sea-surface. With known varying effective dielectric constant of closely packed 
air-bubbles in these layers, estimates of millimetre wave attenuation through layers of sea-foams as 
functions of frequency, foam layer thickness, polarization and angle of incidence are reported. 

Keywords —Extinction, Sea-foams, Split-step Fourier method. 
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I. INTRODUCTION 

Due to the incessant demand for precision and 
accurate prediction of climate and weather variables 
from passive radiometric measurements at WindSat 
frequencies, there is need investigate the direct and 
indirect effects of air-bubbles (sea-foams) 
interactions with millimeter wave at the air- 
seawater interface. Evaluation of extinction of 
millimeter wave by a flat sea-surface covered by 
foam is a significant research problem which 
contributes vastly in satellite-based geophysical 
retrievals of environmental variables such as sea 
surface temperature, sea surface emissivity and 
brightness temperature. 

Anguelova et.al asserted that bubble plumes in 
seawater forms foam layers at the sea surface and in 
some instances sea spray suspended closely above 
sea surface. The skin depth of sea foam at 
microwave frequency narrows the remote sensing 


investigations of the sea foam to its surface 
expression [1]. Many air-sea interaction processes 
are quantified in terms of whitecap fraction W 
because oceanic whitecaps are the most visible and 
direct way of observing breaking of wind waves in 
the open ocean. Enhanced by breaking waves, 
surface fluxes of momentum, heat, and mass are 
critical for ocean atmosphere coupling and thus 
affect the accuracy of models used to forecast 
weather and predict storm intensification and 
climate change. Whitecap fraction is defined as the 
fraction of a unit sea surface covered by foam [2]. It 
has been traditionally measured by extracting the 
high intensity pixels marking white water in still 
photographs or video images collected from towers, 
ships, and aircrafts. Satellite-based passive remote 
sensing of whitecap fraction is a recent 
development that allows long term, consistent 
observations of white capping on a global scale. 
The remote sensing method relies on change of 
ocean surface emissivity at microwave frequencies 
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(e.g., 6 to 37GHz) due to the presence of sea foam 
on a rough sea surface. These changes at the ocean 
surface are observed from the satellite as brightness 
temperatureT# [3]-[4]. 

The algorithm to obtain W from satellite 
observations of T B was developed at the Naval 
Research Laboratory within the framework of 
WindSat mission. It improved upon the feasibility 
study of this remote sensing technique by using 
independent sources for the input variables of the 
algorithm, physically based models for the 
emissivity of rough sea surface and emissivity of 
foam, improved rain flag, and improved 
atmospheric model necessary for the atmospheric 
correction. The database built with this algorithm 
compiles W for entire year 2006 matched in time 
and space with data for the wind vector, wave field 
(such as significant wave height and peak wave 
period), and environmental parameters (such as sea 
surface temperature and atmospheric stability). This 
data base has proved useful in analysing and 
quantifying the variability of W. Magdalena et.al 
presented an updated algorithm for estimating W 
from WindSat T B data using new sources and 
products for the input variables. This approach 
replaces the originally used QuikSCAT data for 
ocean wind vector are replaced with new wind 
vector fields [3]-[4]. 

Xiaobin Tin et.al found large discrepancies when 
comparing measurements and model simulations as 
wind speed (WS) rise above 12ms" 1 . Over the open 
ocean and for moderately wind speeds (WS s ), the 
Soil Moisture and Ocean Salinity (SMOS) 
brightness temperature ( T B ) was initially consistent 
with T b computations made by theoretical 
prelaunch models implemented in European Space 
Level 2 Ocean Salinity processor [5]. In [5] a new 
approach was proposed using new set of parameters 
for sea wave spectrum and foam coverage model 
that can be used for simulating L-band radiometer 
data over a large range of WS based on the 
deduction of wind induced components from the 
SMOS data. 
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[6]-[7] used Monte Carlo Simulation and the dense 
radiative transfer theory (DMRT) with quasi¬ 
crystalline approximation (QCA) to evaluate direct 
and inverse scattering problems. The Monte Carlo 
simulation use an exact numerical formulation 
based on Foldy-Lax multiple scattering equations of 
3-D Maxwell’s equation. Previous calculations 
assumed particles with size parameter (ka « 1) 
using DMRT. An extension to moderate sized 
particle is vital for millimetre wave remote sensing 
because high frequency range above 10GHz, the 
particle size in geographical media are comparable 
to the wavelength (Mie Scattering). Contributions 
by Chi-Te Chan et.al shows that analytical results 
were consistent with Monte Carlo simulations of 
exact solutions of Maxwell’s equations of randomly 
distributed finite size sphere without adjusting 
parameters. Input parameters of the model are all 
physical parameters of sizes, concentrations 
(volume fractions) and permittivity. The results of 
the models were used to obtain brightness 
temperatures in passive remote sensing of ocean 
foam at 19GHz and 37GHz. 

Mie scattering could be used with QCA and 
QCA-CP (quasi-crystalline approximation with 
coherent potential) to provide solution that is 
consistent for particle sizes comparable to or larger 
than the wavelength. With results obtained from [7], 
the extinction behaviour of sea foam was illustrated, 
thermal emission from sea foam was evaluated and 
it was shown that the extinction was dominated by 
absorption. A physical model of foam emission was 
obtained that relates observed brightness 
temperature to the microstructure of foam as well as 
ocean surface wind vector. The brightness 
temperature of sea form can be presented as a 
function of observed angle and frequency and 
radiative transfer equation could be solved using 
derived QCA parameters. 

II. ESTIMATE OF THE EFFECTIVE 
DIELECTRIC CONSTANT OF SEA- 
FORMS 

Figures A-E of [8] comprise of 3D (three 
dimensions) slices of randomly packed spheres in a 
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unit cube. These slices were later translated to 2D 
slices of randomly packed circles. The conversion 
of 3D randomly packed spheres to 2D packed circle 
was achieved by calculating the radii of each 
individual circle using the concept of intersection of 
a sphere and a plane. The 2D slices were discretized 
with grid sizes Ax and Ay which leads to 
intersection of the circles bounded in a unit square 
with some grid points. The grid sizes were sampled 
such that the edges of the circles circumference 
which intersects with grid points farther from the 
inner grids bounded by the circles are negligible. 

Estimate of the effective dielectric constant of sea 
foams were made by modelling the randomly 
packed bubbles as concentric circles in 2-D where 
the outer circle is a mixture of air and seawater 
while the inner circle contains about (80 — 95) % 
air, with these estimates we were able to calculate 
the area of the annulus (ring) as the radii of the 
outer circles are known. 

Table 1 and 2 of [8] illustrate dielectric constant 
of seawater at fixed salinity 34 psu and sea surface 
temperature 20°C. The dielectric constant of air is 
taken as 1.00005 + O.OOOOi . The area of the 
circles in each slice was calculated using the total 
number of grid points. The sea-foam effective 
dielectric constants at WindSat frequencies 6.8 GHz , 
18.7 GHz , and 23.8 GHz were calculated for five 
slices of randomly packed spherical air-bubbles that 
are coated with thin layers of seawater. 

TABLE I 

Results for Dielectric constant of sea foam at frequencies of 6.8 GHz, 

18.7 GHz and 23.8 GHz for 5 2-D Slices of randomly packed air-bubbles 
covered with thin-layer of seawater. 


FREQUENCY 


Slice 1 


Slice 2 


Slice 3 


Slice 4 


Slice 5 


6.8 GHz 


1.1294 + 0.1297 i 


1.2032 + 0.1562 i 


1.2122 + 0.153 7i 


1.3060 + 0.2148* 


1.3457 + 0.2357; 


18.7 GHz 


1.0657 + 0.0896; 


1.0937 + 0.1460; 


1.1252 + 0.1652; 


1.2094 + 0.1804; 


1.2547 + 0.1993; 


23.8 GHz 


1.0008 + 0.126 7i 


1.0117 + 0.1526; 


1.0263 + 0.1833; 


1.0427 + 0.2098 


1.0567 + 0.2302; 


III. PARABOLIC WAVE EQUATION 
METHOD 


Parabolic wave equation method is one of the 
numerous computational and numerical 
approximations used in modelling the interaction of 
EM fields with physical objects and the 
environment. Parabolic equation method (PEM) 
provides efficient approximation to the Helmholtz 
wave equation which is derived by decoupling of 
Maxwell’s equations to evaluate radio-wave 
propagation in random and inhomogeneous media. 

Leontovich and Fock [9,10] introduced the 
parabolic wave equation which is very powerful in 
its application for analysing the problem of radio¬ 
wave diffraction around the earth or in the 
atmosphere. [11] Malyuzhinets, Fock and 
Wainstein generalized the PEM based on studies of 
ray-coordinates and transverse diffusion. 
Malyuzhinets unified geometric optics with the 
parabolic approximate method to establish an 
efficient theory of diffraction of obstacles. Tappert 
et.al [9-11] developed the very powerful split-step 
Fourier solution of the parabolic equation for 
solving underwater acoustic problems. 


PEM applications spans from radio-waves 
propagation problems in VHF to millimetre wave 
affected by atmospheric refraction, diffraction and 
reflection by irregular terrains and rough sea 
surface. More reviews and references of the 
parabolic equation method can be accessed in [11] 
and [9-12]. 

The principal advantages of the various parabolic 
wave equations derived below is that it constitutes 
an initial value problem in range and hence can be 
solved by a range matching numerical technique, 
given a source field distribution over depth at the 
initial range. 

Over the years, several different solution techniques 
have been implemented in computer codes [9-12], 
but only the split-step Fourier technique and various 
linite-difference / finite-element techniques have 
gained widespread use in the underwater acoustic 
community. 

Before going into details on the numerical solution 
schemes, let us briefly point out some advantages 
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and disadvantages of these two main solution 
techniques. 

The split-step algorithm has been extensively used 
to solve the SPE ever since it was developed by 
Hardin and Tappert [9] in the early 1970s. The 
technique is computationally efficient for long- 
range, narrow-angle propagation problems with 
negligible bottom interactions. 

For short-range, deep water problems and shallow 
water problems in general, propagation is basically 
more wide-angled and bottom interacting paths 
become more important. 

This requires the use of wide-angle PEs, which can 
be solved only by finite difference or finite 
elements. Moreover, the strong speed and density 
contrasts encountered at the water-bottom interface 
adversely affect the computational efficiency of the 
split-step technique, which in cases of strong 
bottom interaction requires an excessively fine 
computational grid(Ax,z). Hence, the advantage of 
higher computational efficiency of the split-step 
technique is entirely lost in situations with strong 
bottom interactions. 

Finite-difference and Finite-element solutions are 
applicable to PE for arbitrary large angles. The 
main drawback of these schemes is that, for long 
range the split-step solution is more efficient and 
also for narrow angle with minimal or no bottom 
interaction. The split-step solution remains the most 
adopted technique for performance prediction as it 
is more suitable to solving many practical ocean 
surface problems. 

Conversely, the finite-difference and finite-element 
schemes have widespread application for wide- 
angle and bottom interacting boundaries. It is 
prominent for providing higher accuracy in these 
domains. 

The most recent development in terms of efficient 
PE solution schemes is a split-step Pade’ 
approximations derived by Collins [13,14]. He uses 
higher order Pade’ approximations not the square 
root operators. The result is a considerable 
efficiency gain through the use of higher range step. 
Thus, the scheme is claimed to be more than an 
order of magnitude faster than standard FD/FE 
solution techniques. This could create a unified PE 
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solution approach where the accurate high angle 
PEs can be solved with the efficiency of the 
classical split-step Fourier scheme. 


IV. DERIVATION OF THE PARABOLIC 
WAVE EQUATION 

Using a simple model which describe the 
propagation of a reduced function 
i|;(x,z) = u(x f z)e lkx (\) 
associated with the direction of propagation x. 
Where u(x, z) can be expressed as 
u(x,z) = X ¥(x,z)e~ lkx ( 2) 

The Helmholtz equation of the reduced function is 
obtained by decoupling Maxwell’s equations and 
can be expressed as 
V 2 ¥ + fc 2 n 2 (x,z)¥ = 0(3) 

Where n 2 (x, z) is the refractive index and k is the 
wavenumber. The refractive index n 2 (x, z) is 
assumed to possess smooth variations. The reduced 
function implies that the propagation energy varies 
slowly at angles close to the paraxial direction. 

The Laplacian of ^(x, z) can be expressed as 
V 2l P = [V 2 u + 2 ikVu - k 2 u\e ikx ( 4) 

Where V is taken as the product of a plane wave 
solution, substituting equation (4) into equation (3) 
yields 

[ V 2 u + 2 ikVu — k 2 u]e lkx + 
k 2 n 2 {x,z ) ue lkx = 0(5) 

V 2 u + 2 ikVu + k 2 n 2 {(x,z) — 1 }u = 0(6) 

The Laplacian operator for 2D can be expressed as 

a 2 a 2 

(— + —), in the atmosphere where n — 1 is small, 

a 2 

we neglect — as small (paraxial approximation) 
and equation (6) becomes 

+ 2 ik^ + k 2 (n 2 — 1 )u = 0 (7) 

dz A dx v y 

Where equation (7) is the standard parabolic 
equation. 
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Equation (7) by method of separation of variables 
can be written as 

2ik= (1 - n 2 )k 2 u - j-j(8) 

OX oz z 

Which can be expressed as 

i = '{d -” 2) ‘ 2 -S"( 9) 

V 2 

Let consider A — (1 — n 2 )k 2 -and substitute 

v ' 2 k 

A into equation (9). This yields 
P = iAu( 10) 

dx 

Taking like terms 
p = iAdx( 11) 

The solution of the parabolic wave equation 
becomes 

u(x 0 + A x,z) = u(x 0 )e lAAx ( 12) 

It is worthy of note that equation (7) can be factored 
out to obtain 

{£+ik(l-&}{£;+W + Q)}(13) 

This gives us 

f x = ~ «< 14 > 

y x = -ik(l + 0(15) 

Where Q is the pseudo-differential operator and is 
defined by 

Q = Jp^ + n2 (*' z ) (16) 

Equation (14) is the outgoing parabolic equation 
and equation (15) is the incoming parabolic wave 
equation. 

Considering the propagation medium as 
homogeneous with refractive index n, the field 
component M J satisfies two-dimensional scalar wave 
equation 

Jz* + + k 2 n 2 (z,x) = 0(16) 

The refractive index varies with the range x and 
height z and equation (16) is a good approximation 
provided n varies slowly with wavelength. It is 
worth noting that equation (16) is not exact [Levy]. 
If the propagation medium is vacuum, the standard 
parabolic wave equation in (7) 

24+2^=0(17) 

dz z dx 

The solution of equation (1.43) can be expressed as 

■p 2 . 

u(x 0 + A x,z) = u(x 0 )e‘ 2 k (18) 
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V. THE SPLIT STEP FOURIER TRANSFORM 
SOLUTION 

The split-step Fourier method is a very efficient 
PEM which separate the refractive effect from the 
diffractive part of the propagator. Considering a 
two-dimensional scalar wave equation for 
horizontally and vertically polarised wave. Hardin 
et.al introduced the split-step Fourier method which 
transforms the rough surface problem with 
propagation through a sequence of phase screens 
[14-17], 

The standard parabolic equation (SPE) in equation 
(7) can be written as 

5V = - !)} “09) 


Let 

^ 


B = n 2 {x,z) - 1(21) 


Equation (19) becomes 
d f x = ^{A + B}u( 22) 

The analytic solution of the SPE is 

lie 

u(x + Ax,z) = u(x,z)eT Ax ( A+ ^ (23) 

Using 

5 = ^(24) 

Equation (24) yields 

u(x + Ax,z)= u(x,z)e s ( A+B X 25) 

Equation (25) is the split-step solution which 
represent the field propagating through series of 
phase screens. The field is first propagated through 
a slice of homogeneous medium characterised by 
the exponent of A. 

The SSFM routine was implemented to propagate 
the plane wave 

E(z 0 ,x,y) = E(z,x,y)exp{ik x + ik y + ik z z ) ( 26 ) 
with E(z, x, y) ~ 1 along the forward +z direction. 
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The plane wave was propagated through five (5) 2D 
slices of sea foam layers each containing 
isotropically distributed bubbles. The slices are 
equally dimensioned with area 100mm x 100mm 
with layer thickness 5^ = 10 mm separating adjacent 

layers. The foam layer thickness d » Xq is required 

to account for attenuation (E-field amplitude 
variation) and diffuse scattering (E-field phase 
variation) as the incident E-field travels through 
slices of the sea foam layer. WindSat frequency 
channels (6.8 GHz, 18.7 GHz, 23.8 GHz) were used 
for propagation of the E-field through slices of sea 
foam layers. 



Fig. 1 Geometric configuration of foam-covered sea-surface divided into 
slices of sea-foam layers. 


The incident wave is tilted from the normal so that 
there is an initial phase gradient along the surface of 
the sea foam model. This is done by assigning 

values of 0 < 0j< 90°for the zenith angle and fixed 
azimuthal angle for example cpp 0°or 45°. 

For each of the angles 0jwe compute the E-fields 

that emerges from the foam layer and calculate the 
FFTs of the fields. The contour plots of the 
intensity of the scattered field is displayed. Range 
of the depths of foam for realistic scenarios were 
observed to investigate the influence of the foam as 
a thin phase scattering screen and deep phase 
scattering screen. 
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VI. MILLIMETER WAVE ATTENUATION 

AT LOW AND WINDSAT FREQUENCIES 
AS A FUNCTION OF DEPTH OF SEA- 
FORM 

The attenuation of the field intensity in decibel 
(dB) at a low frequency of 6.8 GHz and sea foam 
slice thickness S t = 0.1 mm which describes a thin 
phase scattering screen, as phase perturbation 
increases as the depth of sea foam layer increases. 
Fig. 2 and 3 illustrate the variation of attenuation of 
the E-field in dB with depth in mm of sea-foam 
layer. This is due to increase in diffuse reflections 
or scattering as the EM wave interacts with 
randomly oriented sea-foams with different 
refractive indices. This happens for both TE and 
TM polarized E-fields. 



Fig. 2 Attenuation of Field Intensity for horizontal polarization (TE) with 
zenith 9i = 45° and azimuth 0 = 0° with depth of sea-foam at 6.8 GHz and 
foam layer thickness S t = 0.1 mm. 



Depth of Seafoam in mm 

Fig. 3 Attenuation of Field Intensity for vertical polarization (TM) with 
zenith = 45° and azimuth 0 = 0° with depth of sea-foam at 6.8 GHz and 
foam layer thickness 8 t = 0.1 mm. 
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Plot of Attenuation in dB against Depth in mm for TE POL 



Fig. 4 Attenuation of Field Intensity for horizontal polarization (TE) with 
zenith Q t = 45° and azimuth 0 = 0° with depth of sea-foam at 18.7 GHz and 
foam layer thickness 8 t = 0.1 mm. 

Similarly, the attenuation of the field intensity in 
decibel (dB) at WindSat frequency of 18.7 GHz and 
sea-foam slice thickness S t — 0.1 mm increases as 
the depth of sea foam layer increases. Fig. 4 and 5 
illustrate the variation of attenuation of the E-field 
in dB with depth in mm of sea-foam layer. This is 
also due to increase in diffuse reflections or 
scattering as the EM wave interacts with randomly 
oriented sea-foams with different refractive indices. 
This happens for both TE and TM polarized E- 
fields. 



Fig. 6 Attenuation of Field Intensity for horizontal polarization (TE) with 
zenith = 45° and azimuth 0 = 0° with depth of sea-foam at 23.8 GHz and 
foam layer thickness 8 t = 0.1 mm. 


Plot of Attenuation in dB against Depth in mm TM POL 



foam layer thickness 8 t = 0.1 mm. 


Plot of Attenuation in dB against Depth in mm TM POL 



foam layer thickness 8 t = 0.1 mm. 


For thin phase scattering screens S t = 0.1 mm 
with a given zenith angle 6 t = 45° and azimuth 
<p = 0°, the attenuation cr(dB) increases with depth 
(mm) of sea-foam and frequency at 23.8 GHz. The 
attenuation at 18.7 GHz is less than that at 
23.8 GHz , which agrees with the fact that there is 
more interaction of the E-field with the particles at 
high frequency. In Fig. 6 and 7, attenuation of the 
E-fields for both TE and TM polarizations are due 
to phase perturbations due to diffuse scattering of 
disordered air-bubbles with varying dielectric 
constants 


VII. CONCLUSIONS 

We have shown from our results that attenuation 
of the E-field is due to diffuse scattering when the 
incident E-field is multiple reflected by sea-foam as 
it propagates through the various slices of sea-foam 
layer. The E-field travels through different paths 
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and is scattered due to re-radiation of the dipole 
moments of the disoriented scatterers which 
absorbs some of the incident E-field and re-radiates 
it at various directions. The varying effective 
dielectric constants of the particles enhances diffuse 
scattering of the incident and transmitted E-fields 
within the sea-foam layer. The split-step Fourier 
transform method is a range marching technique 
that is suitable for evaluation of the E-field 
propagation given an initial range, the expected 
depth of the field can be computed as a function of 
range and depth. Future works will evaluate 
variation of the field intensity due to its interaction 
with rough sea-surface covered by foams. 
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